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FOR  GENERAL  NON-LINEAR  SYSTEMS  (U) 


by 
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SUMMARY 

A  computer  program  has  been  developed  for  the  Maximum  Likelihood 
estimation  of  parameters  in  general  non-linear  systems.  Sensitivity  matrix 
elements  are  calculated  numerically,  overcoming  the  need  for  explicit 
sensitivity  equations.  Parameters  such  as  break  points  and  time  shifts  are 
successfully  determined  using  both  simulated  and  actual  test  data. 


DSTO^ 

M  E  I  B  0  U  R  n'e 


(O  COMMONWEALTH  OF  AUSTRALIA  1988 

POSTAL  ADDRESS:  Director,  Aeronautical  Research  Laboratory, 

P.O.  Box  4331,  Melbourne,  Victoria,  3001,  Australia 


CONTENTS 


Page  No. 

'NOTATION  . ii 

1.  INTRODUCTION . 1 

2.  DESCEUPnON  OF  METHOD . . . 1 

3.  EXAMPLES  USING  SIMULATED  DATA . 3 

3.1  Validation  of  Numerical  Sensitivity  Matrix  Computation . 3 

3.1.1  Explicit  Sensitivity  Matrix  Calculation . 4 

3.1.2  Numerical  Sensitivity  Matrix  Computation . 5 

3.2  Systems  with  Break  Points . 6 

3.3  Systems  with  Time  Shifts . 8 

4.  RESULTS  USING  FLIGHT  TEST  DATA . 9 

5.  CONCLUDING  REMARKS . 11 

ACKNOWLEDGEMENT. . 12 


REFERENCES  . 

APPENDIX  A  . 

FIGURES 

DISTRIBUTION 

DOCUMENT  CONTROL  DATA 


.13 

.14 


1  Acosssloo  For  j 

RTIS 

DTIC 

Justl 

6RAAI  df 
rAB  O 
suDoed  □ 

Distrlhutton/ 

Availability  Codes 

Dist 

Ai 

Avail 

Spec 

nad/or 

lal 

. . . 


NOTATION 


GENERAL ; 


dt 

r% 

Time  interval  between  succes^ve  time  measurements,  tj^t  *  ti 

s/ 

n 

Measurement  noise  vector 

N 

Number  of  discrete  time  points 

t 

Time 

f 

Time  prior  to  time  shifting 

ti 

i’*'  discrete  time  point 

u 

Control  input  vector 

X 

State  vector 

y 

Observation  vector 

yi 

jti'  observation  variable 

z 

Measurement  vector 

8^ 

Increment  in  ^  used  for  numerical  sensitivity  calculations 

Change  in  \  per  Keration 

% 

Parameter  vector 

k**  parameter 

v^y 

Sensitivity  matrix 

^^yjk 

jk**  element  of  sensitivity  matrix  (-dy/d^ 

LANDING  GEAR  DROP  TESTS  (SECTION  3) 

Ci 

Tyre  'spring*  constant 

d 

Oleo  deflection 

^mw 

Maximum  value  of  d 

d* 

Tyre  compression 

do 

Break  point 

dt,  da 

Rrst,  second  stage  oleo  deflection  (two  stage  type) 

g 

Acceleration  due  to  gravity 

Gt.Ga 

Oleo  damping  terms 

Kt.Ka 

Oleo  "spring”  constants 

L 

Load  on  landing  gear 

Lo.L. 

Load  on  oleo,  tyre 

M 

Landing  gear  mass 

w 

d+cT 

Wo 

initial  drop  velocity 

Td.tt 

Time  shifts  in  d.  L 

FLIGHT  TEST  (SECTION  4) 

an 

Normal  acceleration 

c 

Reference  chord 

Cm 

Moment  coefficient,  as  a  function  of  Cm^  Cm^,  Cm,.  Cm^.  and  Cm„ 

Cn 

Force  coefficient,  as  a  function  of  Cn^,  C|^,  Cn,,.  and  Cn^ 

i 


g  Acceleration  due  to  gravity 

ly  Moment  of  inertia  sdsout  pitch  axis 
m  Aircraft  mass 

q  Rate  of  change  of  angle  of  pitch 
q  Dynamic  pressure 

R  Radian  to  degree  conversion  factor 

S  Wing  area 

V  Aircraft  velocity 

Xq,  Xa„  Longitudinal  instrument  offsets  from  c.  g. 
Zap  Vertical  instrument  offset  from  c.  g. 

a  Angle  of  attack 

5  Elevator  deflection 

d  Angle  of  pitch 


1.  INTRODUCTION 


The  Maximum  Likelihood  Parameter  Estimation  technique  is  widel/  used  to  determine 
aircrah  flight  parameters  from  flight  test  data  (Ref.  1).  Much  effort  in  this  field  has  been 
confined  to  analysis  of  linear  systems.  If  the  model  is  non-linear,  the  problem  becomes  more 
difficult. 

Here  at  Aeronautical  Research  Laboratory  (ARL),  a  Maximum  Likelihood  computer  program 
was  developed  (Ref.  2)  to  solve  such  non-linear  problems.  This  program  was  successfully 
used  to  determine  accelerometer  offsets  and  calibration  errors,  given  dynamic  flight  test  data 
(Ref.  3).  However,  a  sizable  proportion  of  the  program  was  problem-specific,  in  particular  the 
evaluation  of  the  sensitivity  matrix.  Sensitivity  matrix  elements  were  calculated  explicitly  by 
mathematical  differentiation  of  the  state  equations,  often  a  long  and  tedious  process. 

In  this  document,  a  Maximum  Likelihood  Parameter  Estimation  computer  program  for 
general  non-linear  systems  is  described.  Sensitivities  are  calculated  numerically  by  finite 
differences,  overcoming  the  need  for  explicit  sensitivity  equations.  This  method  allows  the 
user  to  select,  as  unknown  parameters,  quantities  such  as  break  points  and  time  shifts,  for 
which  sensitiv'ities  are  not  known  in  explicit  form. 

Other  Maximum  Likelihood  programs  representing  generalized  non-linear  systems  are 
available  (Ref.  4).  However,  their  specific  application  to  the  estimation  of  break  points  or  time 
shifts  does  not  appear  to  have  been  repotted.  The  program  described  here  is  successfully 
used  to  determine  both  of  these  quantities. 

Section  2  provides  a  brief  theoretical  description  of  the  Maximum  Likelihood  method  for 
non-linear  systems,  and  the  procedure  used  here  to  obtain  sensitivities  numerically.  In 
Section  3,  the  computer  program  developed  is  validated  using  simulated  data,  before  being 
applied  to  a  study  of  aircraft  landing  gear  modelling,  a  topic  currently  being  investigated  by  the 
Aircraft  Behaviour  Studies  -  Rotary  Wing  Group  at  ARL.  In  Section  4,  the  program  is  used  to 
estimate  time  shifts  in  measured  flight  test  data. 


2.  DESCRIPTION  OF  METHOD 


Assume  that  the  system  can  be  described  in  general  by  a  set  of  non-linear  dynamic 
equations  of  the  form ; 

i(t)  -  f{  x(t).  u(t),  o  m 

y(t)  -  g(x(t),u(t),u  (21 

z(ti)  -  y(t()  +  n(t|)  (3] 


where : 

X  is  the  state  vector 
u  is  the  control  input  vector 


■  y  is  the  observation  vector 


z  is  the  measurement  vector,  sampled  at  N  discrete  time  points,  ti ,  for  i  - 1 N 
n  is  the  measurement  noise  vector,  assumed  to  be  Gaussian  with  zero  mean 
4  ^  the  vector  of  unknown  parameters 

The  Maximum  Likelihood  method  (described  in  Ref.  2)  determines  the  most  probable  value 
of  ^  by  an  iterative  procedure  which  can  be  summarized  as  follows : 

R  -  -  y(ti)F 

ia1 

N  N 

i=1  i«1 

where  R  is  the  covariance  of  residuals  and  is  the  change  in  ^  per  iteration. 

Given  R  and  we  obtain  and  hence  an  improved  value  of  which  is  used  to  obtain  a 
new  y(ti)  and  thus  improved  R.  This  process  is  repeated  untii  convergence  is  achieved. 
Computation  of  requires  at  each  time  point,  t) : 


i)  values  of  the  measurement  vector,  z(ti) 
iO  values  of  the  obsenration  vector,  y(tj) 
ifi)  the  sensitivity  matrix,  V^y(ti) 


Values  of  the  measurement  vector  are  read  in  as  data.  Computation  of  the  current 
observation  vector,  y(ti),  from  [2]  requires  current  state  vector  values,  which  are  obtained  by 
numerical  integration  of  the  assumed  system  state  equations  ( [1]  above  ).  A  fourth-order 
Runge-Kutta  numerical  integration  procedure  is  adopted.  The  sensitivity  matrix  elements  are 
here  approximated  by  numerical  differences.  This  approach  overcomes  the  need  for  explicit 
sensitivity  equations,  which  are  not  always  easy  to  determine.  The  central  difference  method 
is  adopted  here,  requiring  evaluation  of  state  and  observation  variables  at  two  perturbed 
parameter  values,  ^+5^  and  %-S^.  The  jk*  element  of  the  sensitivity  matrix,  v^yltOj^,  is  given 
by: 


v^ylWjK 


axK  "  25^ 


[61 


where  yj  and  ^  are  components  of  vectors  y  and 

The  program  developed  here  is  based  on  an  earlier  ARL  Maximum  Likelihood  program  (Ref. 
2).  Details  of  program  changes  are  listed  in  Appendix  A. 


3.  EXAMPLES  USING  SIMULATED  DATA 


Simulated  data  is  used  to  examine  a  number  of  different  systems.  In  3.1 ,  the  numerical 
sensitivity  matrix  computation  is  validated  by  analysing  a  simple  system,  both  with  explicit  and 
numerical  sensitivity  matrix  computations.  In  3.2  and  3.3,  the  numerical  model  is  used  to 
examine  more  complicated  systems  which  include  the  effect  of  break  points  and  time  shifts. 


3.1  Validation  of  Numerical  Sensitivity  Matrix  Computation 

Consider  the  case  of  an  aircraft  landing  gear  drop  test.  The  landing  gear,  comprising  a  large 
mass,  M,  attached  to  an  oleo  and  tyre,  is  dropped  with  initial  velocity  wq.  The  oleo  is  modelled 
as  a  massless,  non-linear  damped  spring,  and  the  tyre  as  a  massless,  linear  undamped  spring. 

The  load  on  the  oleo,  Lo,  is  given  by : 

U  »  K,d2  +  G,d  [7] 

where  d  is  the  oleo  deflection.  The  load  on  the  tyre,  Lt  ,  is  given  by,: 

L,  »  C,d*  [8] 

where  d*  is  the  tyre  compression.  Since  the  tyre  and  oleo  are  assumed  to  be  massless,  the 
oleo  load  is  equal  to  the  tyre  load,  thus  Lo  =  L  ^  I-  (see  Fig.  1). 


The  equation  of  motion  for  the  system  is : 

M(d+d*)  »  Mg-L  [9] 

Taking  d,  d*,  and  w  (=d  +  d*)  as  state  variables,  the  following  state  equations  are  obtained : 

[10] 

Cid*  -  Kid2 


Cid* 

w  =  g  — 


■ST 


d*  =  w-d 


[111 

1121 


Observation  variables  are  the  load  and  oleo  deflection,  given  by  : 

Lotw  =■  Cid*  [13] 

dob*  *  d  [14] 

There  are  no  control  inputs,  u,  in  this  example.  The  unknown  parameters  are  Ki,  Ci.  and 
Gi- 


The  above  equations,  [10]  -  [14] ,  are  subject  to  constraints  that  d,  d*,  and  L  are  all  ^  0. 

The  landing  gear  is  dropped  with  initial  velocity  Wo ,  so  initial  values  of  w,  d,  and  d*  are  Wq,  0, 
and  0  respectively. 


*  Using  simulated  time  histories  of  oieo  toad  and  deflection,  with  81  data  points  and 
At  0.01  s,  the  Maximum  Likelihood  method  is  applied  with  a)  explicit  sensitivity  matrix 
equations  and  b)  numerical  sensitivity  matrix  computations.  Zero  mean  noise  with  an  RMS  of 


0.0025  m  for  d  and  0.5  kN  for  L  is  superimposed  on  the  simulated  data. 

0 

3.1.1  Explicit  Sensitivity  Matrix  Caicuiation 

Sensitivity  matrix  elements,  3yi/3§k.  are  given  by  dki,  dci,  dci,  Lk,  (»  CidKi*),  Lci 
(*Cidci*+d*)  and  Lqi  (=  Cidc,*)  where  dKi  =»  3d/3Ki  etc.  Sensitivities  are  obtained  from  partial 

derivatives  of  the  state  equations  [10]  -  [12] : 

dKi  =  ^(C,  dKi*  -  2dKidKi  -d2) 

[15] 

[16] 

Wki  « 

[17] 

dc,  =  ^(d*  +  C,dci*  -  2dKidci  ) 

[18] 

dci*  -  -  dci 

[19] 

wci  =  ^(C,dc,*+d2) 

[20] 

d<3,  -  ^(^C,d<3,*  -  2dK,dQ,  - 

)  [21] 

dci*  »  Wsi  -  dg, 

[22] 

-Cldc,* 

Wg,  -  ^ 

[23] 

The  above  derivatives  are  numerically  integrated  with  respect  to  time  using  a  fourth-order 
Runge  Kutta  procedure,  and  the  sensitivity  elements  calculated,  from  which  the  most 
probable  values  for  parameters  K,,  C,,  and  Gi  are  determined.  ^ 

After  10  iterations,  excellent  agreement  with  the  simulated  data  is  obtained  (Fig.  2).  The 
parameter  values  are  listed  in  Table  1  with  the  true  values  being  the  values  used  in  the 
simulation,  and  the  a  priori  values  being  the  initial  guess  of  these  values.  The  Cramer  -  Rao 
error  bounds  are  also  shown. 


Table  1.  Parameter  Estimates  for  3-Parameter  Model 


Parameter 

Unit 

A  Priori 

Value 

Taie 

Value 

Maximum  Likelihood  Parameter  Values 

Explicit  Sensitivities 

Numerical  Sensitivities 

K, 

Nfrr2(x10®) 

1 

4 

3.993±0.009 

3.993±0.009 

Gi 

Nsrrr’  (xIO*) 

1 

2.5 

2.496±0.007 

2.496±0.007 

c, 

Nntr’  (xIQS) 

1 

7 

7.037±0.041 

7.038±0.041 

RunTime  (s) 

4.29 

11.62 

3.1.2  Numerical  Sensitivity  Matrix  Computation 

»  » 

Sensitivity  matrix  elements  are  calculated  using  the  central  difference  method  {see  [6] )  with 
=  10-3-  After  10  iterations,  the  Maximum  Likelihood  method  results  in  the  parameter 
values  shown  in  Table  1.  It  is  seen  that  results  are  almost  identical  to  the  explicit  sensitivity 
case,  and  indeed  graphs  of  oleo  deflection  and  load  versus  time  (Rg.  3)  are  found  to  be 
indistinguishable  from  the  expBcit  sensitivity  graphs  (Fig.  2).  There  is  however  an  increase  in 
computer  mn  time  due  to  the  numerical  integration  of  the  state  equations  for  incremented 
parameters  §  +  5?  and  §  -  54.  as  well  as  the  usual  4* 

Other  test  cases  (not  shown  here)  were  also  found  to  give  excellent  agreement  between 
explicit  and  numerical  sensitivity  calculations,  thus  validating  the  numerical  procedure  used 
here  for  the  computation  of  the  sensitivity  elements. 

The  effect  of  varying  the  size  of  54k/4k  in  the  numerical  sensitivity  case  was  examined.  For 
54k/4k  ranging  from  0.1  down  to  10*^,  no  appreciable  increase  in  accuracy  was  obtained  for  the 
simulated  data  of  Rg.  2  (see  Table  2).  >■ 


Table  2,  Maximum  Likelihood  Parameter  Estimates  using  Numerical  Sensitivity 
Calculations  with  Varying 


Parameter 


Unit 

A  Priori 

Tme 

5CK/4k 

Value 

Value 

0.1 

0.01 

0.001 

Nnr2(x10®) 

1 

4 

3.9928 

3.9928 

3.9928 

Nsrrr’  (xIO*) 

1 

2.5 

2.4962 

2.4962 

2.4962 

Nnr’  (x10®) 

1 

7 

7.0380 

7.0376 

7.0376 

3.2  Systems  with  Break  Points 

A  more  acfvanced  landing  gear  system  is  the  two-stage  type.  For  oleo  deflection,  d,  greater 
than  some  value,  do,  a  second,  stiffer,  more  damped  'spring*  is  activated.  In  this  way,  extra 
hard  landings  are  catered  for,  whilst  soft  landings  do  not  suffer  fK)m  overdamping.  Essentially, 
a  different  set  of  state  equations  exist  for  d  >  do  with  do  termed  as  a  break  point. 


Assume  the  oleo  load  to  be  modelled  as 

The  state  equations  in  this  case  are : 

Cid* 

w  -  g 

Cid*  -  K,di2- KjdzZ 
d - g - 

d*  -  w-d 


K,d2  +  Gid 

Kido2  +  K2(d-do)2  +  Gad 


(for  d  <  do) 
(for  d  a  do) 


where : 


d<ct 


do  '] 

d-do  ^  d  s  do 

Ga  J 


di  is  the  first  stage  deflection  and  da  the  second  stage  deflection.  Observation  variables  are 
the  oleo  load  and  deflection. 

Along  with  parameters  Ki,  Ka,  Gi.Ga  and  Ci,  the  break  point  do  is  taken  as  an  unknown 
parameter.  Using  simulated  time  histories  of  oleo*  load  and  deflection,  the  Maximum 
Likelihood  method  with  numerical  sensitivity  calculations  is  applied.  The  calculated 
observations  give  excellent  agreement  with  the  simulated  measurements  (Fig.  3).  The 
parameter  values  after  1 0  iterations  are  given  in  Table  3. 


Table  3.  Parameter  Estimates  for  6-Parameter  Model  (Including  Break  Point) 


Parameter 

Unit 

1 - 

True 

Value 

Maximum  Likelihood 

Parameter  Value 

Ki 

Nm-2(x10®) 

2 

4 

4.040  ±  0.017 

Kz 

Nm-2  (xIO®) 

10 

45 

45.825  ±  1.348 

Gi 

Nsnr'  (xIO^) 

1.5 

2.5 

2.473  ±  0.009 

Gz 

Nsm-^  (x10^) 

3.5 

4 

3.966  ±  0.030 

Cl 

Nm-’  (xIOS) 

4 

7 

7.064  ±  0.039 

do 

m 

0.1 

0.23 

0.227  ±  0.001 

It  should  be  noted  that  careful  consideration  needs  to  be  given  to  the  selection  of  the 
a  priori  value  of  the  break  point  parameter  do-  Clearly,  if  do  >  dmax  {where  dmax  is  the  maximum 
value  reached  by  quantity  d  in  the  time  interval  under  consideration),  either  initially  or  during 
one  of  the  early  iterations  (when  parameter  values  are  liable  to  oscillate),  then  any  small 
change  in  do  will  have  no  effect  on  the  obsenration  vector,  y.  Consequently,  sensitivities 
ayi/3do  will  all  be  zero  resulting  in  no  improved  estimate  of  do.  A  suitable  initial  value  for  do  can 
usually  be  obtained  by  experimentation. 

For  most  parameters,  the  size  of  used  in  the  numerical  sensitivity  calculations,  should  be 
significantly  less  than  the  size  of  We  use  here  5^/^^  =  10"^.  However,  break  points  or  time 
shifts  are  special  cases,  and  careful  consideration  needs  to  be  given  to  the  size  of  5^.  This  is 
because  measurements,  z,  only  exist  at  discrete  time  points,  ti,  and  observations,  y,  are  only 
calculated  at  these  time  points.  For  break  point  d,  (time  shifts  are  discussed  in  Section  3.3),  in 
order  that  3yj/3do  is  non-zero  for  at  least  one  tj,  any  change,  5do,  in  the  value  of  do  must  be 
large  enough  such  that  the  time  where  the  break  point  acts,  t(do),  moves  across  at  least  one 
time  point.  If  t(do)  is  within  the  time  interval  (tj,  tj*,],  then  for  5do  to  result  in  non-zero 
sensitivities,  we  require: 

|t(do)-t(do-5do)|  >  t(do)-ti 

or  lt(do  +  5do)-t(do)|  >  U,-t(do)  {seeFQ.5) 


The  size  of  5do  necessary  to  satisfy  the  above  condition  can  be  found  by  experimentation. 
In  the  case  reported  here,  5do  »  0.01  was  sufficient. 


3.3  Systems  with  Time  Shifts 


When  recording  actual  test  data,  time  shifts  between  different  measurements  can  occur  as  a 
result  of  instrumentation  lag  (Ref.  3).  In  general,  the  time  shift  for  a  particular  measurement  is 
not  known,  and  can  be  included  in  the  Maximum  Likelihood  procedure  as  an  additional 
parameter.  Time  shifts  are  expected  to  be  small  and  a  priori  values  are  usually  set  to  zero. 

In  a  landing  gear  drop  test  such  as  that  in  Section  3.2,  any  time  shifts  in  measurements  d 
and  L  need  not  be  the  same  and  can  be  represented  by  unknown  parameters  T<j  and 
respectively. 

The  simulated  drop  test  data  from  Section  3.2  is  used,  with  time  shifts  added  in.  Assuming 
the  same  landing  gear  model  as  in  Section  3.2,  we  aim  to  determine  these  time  shifts  using 
the  Maximum  Likelihood  procedure. 

The  state  and  observation  equations  ( Section  3.2 )  are : 


c,d*(t )  -  Kidi2(t)  -  K2d22(ti 

[28] 

d*(f)  =  w(t)-d(ti 

[29] 

Mf) = 

[30] 

dob*(t)  -  d(t) 

[31J 

Lob,(t)  =  c,<r(f) 

[32] 

where  f  is  the  time  without  any  time  shift. 

A  time  shift  in  d  or  L  will  result  in  a  translation  of  the  above  observations  ; 

dob,(t)  =  d{t-td)  [33] 

Lote(t)  =  L(t-T0  [34] 

Application  of  the  Maximum  Likelihood  procedure  gives  excellent  results  in  10  iterations 
(Fig.  6  and  Table  4)  with  both  time  shifts  accurately  predicted.  Run  time  has  increased  from 
24.13  to  30.50  seconds;  however,  this  is  to  be  expected  since  the  number  of  unknown 
parameters  has  increased  from  6  to  8. 


Table  4.  Parameter  Estimates  for  8-Parameter  Mode!  (Includhg  Break  Point  and  Time  Shifts) 


parameter 

Unit 

HI 

Troe 

Value 

Maximum  Likelihood 

Parameter  Value 

K, 

Nm-2(x10S) 

2 

4 

4.031  ±  0.028 

Kz 

Nm-2(x10S) 

20 

45 

47.535  ±  1.977 

Gi 

Nsmr’  (xtO^) 

1.5 

2.5 

2.486  ±  0.018 

Gz 

Nsnr^  (x10^) 

3 

4 

3.931  ±  0.054 

Cl 

Nrrr’  (xtO®) 

5 

7 

7.039  ±0.110 

do 

m 

0.1 

0.23 

0.229  ±  0.002 

•Cd 

s 

0 

0.07 

0.070±0.0008 

-Cl 

s 

0 

0.09  * 

0.090±0.0005 

In  the  Maximum  Likelihood  program  developed  here,  any  parameter  representing  a  time 
shift  is  automatically  rounded  to  a  multiple  of  dt,  the  time  intenral  between  successive 
measurements,  tj.  Observations  are  only  computed  at  time  t-,  for  comparison  with  the 
measurements,  and  thus  any  time  shift  must  of  necessity  be  a  multiple  of  dt.  Likewise,  the 
small  increment  in  time  lag  parameter  required  for  the  numerical  evaluation  of  sensitivities, 
must  also  be  a  multiple  of  dt.  Sensitivities  are  evaluated  here  by  examining  two  outputs  that 
have  time  lags  i  +  dt  and  x  -  dt.  Then : 

3y^  yi(T+dt)  -  yi(t-dt) 

3t  “  2dt 

i 

4.  RESULTS  USING  FLIGHT  TEST  DATA 

Fixed  wing  flight  test  data  was  available  for  an  aircraft  undergoing  a  longitudinal  manoeuvre. 
State  variables  are  taken  as  a,  q,  and  6  with  control  input  5.  Observation  variables,  for  which 
measurements  are  given,  are  a.  q,  and  an  where : 

a  -  angle  of  attack  (deg) 

0  >  pitch  angle  (deg) 

q  -  rate  of  change  of  pitch  angle,  0  (deg  S’’) 
an  >  normal  acceleration  (ft  s*^ 

S  -  elevator  deflection  (deg) 
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’  The  full  state  equations  for  a  longitudinal  manoeuvre  are  given  in  Equation  [55],  Ref.  5.  For 

« 

small  a  and  we  obtain  the  state  equations : 


with : 


Cn 


Cm 


CnjjO  +  Cn^  +  CnjS  +  C,^ 

C»^o  +  Cm,^^  +  CM58  + 


[35] 

[36] 

[37] 

(Force  Coefficient) 
(Moment  Coefficient) 


Observations  (Equation  [56],  Ref.  5)  are  : 


®ob* 

dob* 

S)b« 


“-vrq 


“Cn  +  TTrq  +  Trrr 


mg 


gR  R2g ' 


[38] 

[39] 

[40] 


where : 

R  -  Radian  to  degree  conversion  factor  -  57.2958  deg  rach’ 

S  =  Wing  area  =>  550  ft* 
m  «  Aircraft  mass  <=  1929.88  slug 
V  =  Aircraft  velocity,  assumed  constant,  »  1088.4199  ft  s"’ 
g  »  Acceleration  due  to  gravity  =  32.1 74  ft  s'* 
c  =  Reference  chord  =  8.8  ft 

ly  ■  Moment  of  inertia  about  pitch  axis  »  31 0  91 2.906  slug  ft* 

Xa  =  o  -  vane  Longitudinal  offset  from  c.  g.  =  26.3065  ft 
Ztf,  =  Accelerometer  (normal)  vertical  offset  from  c.  g.  =  1 .3685  ft 
X|^,  ■»  Accelerometer  (normal)  longitudinal  offset  from  c.  g.  *  3.371 5  ft 
Cn<,  »  8.75  rad-’ 

Cnj  »  0.015  deg*’ 


i\  »  Dynamic  pressure,  assumed  constant,  *  1318.026  Ibf  ft* 


with  initial  conditions  Oo  -  2.3169  deg.  do  -  0.0948  deg  s•^  and  Qo  -  2.4568  deg.  Unknown 
parameters  are  Cm^,  Cm^.  Cm,,  Cm^,  Cn„.  and  Cn^  +  Oo .  Parameter  Cm.^  is  also  unknown, 
but  in  the  manoeuvre  reported  here,  it  is  difficult  to  determine  both  Cm^  and  Cm^ 
independently.  Consequently,  the  two  parameters  are  linked  by  a  factor,  obtained  from  their 
a  prioii, values : 

■  0.2836  X  Cm^ 

The  Maximum  Likelihood  procedure  is  applied  without  time  shifts,  with  a  time  interval  of 
dt  -  1/60  sec  and  10  iterations,  and  with  the  control  input  variation  shown  in  Rg.  7.  Results 
are  shown  in  Fig.  8  and  Table  5.  It  is  noted  that  there  are  discrepancies  between  actual 
measurements  and  predicted  observations  which  appear  to  be  due  to  time  lags.  Application 
of  the  Maximum  Likelihood  procedure  with  three  additional  time  shift  parameters,  Xq.  and 
leads  to  marked  improvements  in  the  predicted  observation  variables  (Rg.  8).  The 
Cramer-Rao  error  bounds  are  also  reduced  (Table  5)  indicating  improved  parameter 
identification. 


Tables.  Parameter  Estimates  using  Flight  Test  Data 


Parameter 

Unit 

A  Priori 

Value 

Maximum  Likelihood  Parameter  Values 

Without  Time  Shifts 

With  Time  Shifts 

deg*’ 

•0.050 

-0.036  ±  0.0002 

-0.041  ±  0.00008 

Cm, 

rad-’ 

-30 

-9.92  ±  0.46 

-8.29  ±  0.14 

Cmj 

deg*’ 

-0.035 

-0.019  ±  0.0002 

-0.022  ±  0.00009 

Cmo 

- 

0 

0.102  ±  0.0007 

0.116  ±  0.0003 

Ch. 

deg-’ 

0.09 

0.071  ±  0.0009 

0.071  ±  0.0003 

- 

0 

-0.093  ±  0.003 

-0.092  ±  0.0008 

Cua  +  Olo 

- 

-0.25 

-0.087  ±  0.004 

-0.093  ±  0.001 

s 

0 

- 

0.033  ±  0.001 

■tq 

s 

0 

- 

0.050  ±  0.0007 

1 

s 

0 

1 

1 

0.067  ±  0.0007 

5.  CONCLUDING  REMARKS 

A  general  Maximum  Likelihood  program,  suitable  for  identification  of  parameters  in  non¬ 
linear  systems,  has  been  described.  A  summary  of  the  theoretical  background  has  been 
given,  and  differences  between  an  earlier  ARL  program  noted.  Sensitivities  are  now 
calculated  numerically,  a  process  which  allows  time  shifts  or  break  points  to  be  identified.  The 
program  was  validated  on  a  simple  3-parameter  system,  before  being  applied  to  a  study  of 
aircraft  landing  gear  modelling.  The  study,  using  simulated  data,  has  shown  that  the  program 


can  successfully  identify  time  shifts,  break  points,  and  conventional  parameters,  using  the 
numerical  sensitivity  approach.  Finally,  the  program  has  been  applied  to  actual  flight  test  data, 
and  has  successfully  matched  time  shifts  and  conventional  parameters. 
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APPENDIX  A 


Computer  Program  Structure 

Th^  program  described  in  this  document  was  developed  from  an  earlier  ARL  program, 
COMP  AT  (Ref.  2).  In  its  earlier  form,  COMP  AT  required  four  problem-dependent 
'subroutines,  which  calculated : 

•  Initial  conditions 

•  Explicit  equations  for  derivatives  of  sensitivities 

•  Derivatives  of  state  equations 

•  Output  responses  (observations) 

In  the  modified  program,  COMPAT.JB.5,  initial  conditions  are  now  specified  in  the  datafile, 

COMDAT.JB,  and  sensitivity  matrix  elements  are  evaluated  numerically,  requiring  only  two 
problem-dependent  subroutines  (RESPJB  and  DERIV.JB.3  which  calculate  the  output 
responses  and  state  derivative  equations  respectively).  The  final*  program  stnjcture  is  : 

Main  program:  COMPAT.JB.5 

Subroutines  :  RESP.JB  (output  response  equations) 

DERI\/.JB.3  (derivatives  of  state  vector  equations) 

COM.SUB.JB  (Maximum  Likelihood  iterative  procedure) 

Data  file :  COMDAT.JB  (input  the  necessary  information  to  specify  the 

problem,  eg.  number  of  states,  number  of 
parameters ,  as  well  as  a  priori  parameter  values 
and  initial  state  conditions.  Rnaily,  time  histories 
of  control  inputs,  u(tj),  and  measurements, 
z(ti).  are  read  in) 

Output  files  :  HP1.PLOT  (contains  time  histories  of  the  inputs,  and 

measured  &  calculated  outputs  suitable  for  * 
producing  time  history  plots) 

COMOUT.JB.3  (contains  an  iterative  history,  as  well  as  final 

values  of  parameters  and  their  Cramer-Rao 
error  bounds) 
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Figure  2.  Oleo  Load  and  Displacement  (or  3-Parameter  Model  (Explicit  Sensitivity 
Equations)  w,  ■  4  ms  \  M-2000  kg 
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Figure  3.  Oleo  Load  and  Displacement  for  3-Parameter  Model  (Numerical  Sensitivity 
Computations)  w,  -  4  ms'L  M-2000  kg 
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Figure  4.  Oleo  Load  and  Displacement  for  6-Parameter  Model  (Including  Break  Point ) 

Wo  -  4  msr’ ,  M-2000  kg 


Figure  7.  Control  Input  (Elevator  Deflection)  Variation  with  Time  -  Flight  Test 


Figure  8.  Variation  of  Angle  of  Attack,  Pitch  Rate,  and  Normal  Acceleration  -  Flight  Test 
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